if (!exists("app_root")) {
  stop("Need to have app_root defined to source this template");
}

# hack-hack-hack
library("Biobase");
library("Hmisc");

source ( paste ( app_root, "cgi-lib/R/ep.io.R", sep="/" ) );
source ( paste ( app_root, "cgi-lib/R/ep.statistics.R", sep="/" ) );

ep.ttest <- function ( src, src2, classone, classtwo, alpha, correction, controlgeneids, controlmean, outputfile ) {
  M <- ep.readBin ( src );
  nasM <- apply ( is.na (M), 1, sum );
  Mgood <- M[nasM/ncol(M)<0.20,];

  if ( length(controlgeneids) != 0 ) {
    controlmean <- mean ( M[controlgeneids,], na.rm=T );
    controlsd <- sd(as.vector(M[controlgeneids,]),na.rm=T);
    w <- which ( nasM/ncol(M) < 0.20 );
    w2 <- setdiff ( w, controlgeneids );
    Mgood <- M[w2,];
  } else {
    controlmean <- 0;
    controlsd <- 0;
  }

  if ( src2 != "" ) {
    M2 <- ep.readBin ( src2 );
    nasM2 <- apply ( is.na (M2), 1, sum );
    M2good <- M2[nasM2/ncol(M2)<0.20,];

    if ( nrow(Mgood) != nrow (M2good) ) stop ( "The two datasets need to have the same number of rows with no more than 20% missing values!" );

    pvals <- t ( get.pval.ttest.twosample ( cbind ( Mgood, M2good ), 0:(ncol(Mgood)-1), ncol(Mgood):(ncol(Mgood)+ncol(M2good)-1) ) );
    rownames(pvals) <- rownames(Mgood);
  } else {
    if ( length ( classone ) == 0 && length ( classtwo ) == 0 ) {
      pvals <- t ( get.pval.ttest.onesample (Mgood, Mean = controlmean ) );
      rownames(pvals) <- rownames(Mgood);
    } else {
      pvals <- t ( get.pval.ttest.twosample (Mgood, classone, classtwo ) );
      rownames(pvals) <- rownames(Mgood);
    }
  }

  cutoff <- 15;

  if ( correction == "bonferroni" ) {
    pvals[,1] <- p.adjust ( pvals[,1], "bonferroni" )$adjp;
  } else if ( correction == "hochberg" ) {
    pvals[,1] <- p.adjust ( pvals[,1], "hochberg" )$adjp;
  } else if ( correction == "holm" ) {
    pvals[,1] <- p.adjust ( pvals[,1], "holm" )$adjp;
  }

  orders <- order ( pvals[,1] );
  opvals <- as.data.frame(pvals[orders,]);
  q <- opvals[1:cutoff,];

  png(paste(outputfile, "_errbar.png", sep=""), width=600, pointsize=8);
  errbar_range <- 1.2 * max ( abs ( controlmean ), max ( abs ( q ) ), abs(controlsd) + abs(controlmean) );

  plot.new();
  errbar ( paste(format.pval(q[cutoff:1,1]),' (',rownames(q[cutoff:1,]),')',sep=''), q[cutoff:1,3], q[cutoff:1,2], q[cutoff:1,4], ylim=c(-errbar_range,errbar_range), xlim=c(0,max(max(q[,1]),1.5*alpha) ), ylab=paste('expression level (mean/difference of means)\n(red: p < ',alpha,')',sep='') );

  title ( paste ( cutoff,' lowest p-values with 95% confidence intervals\n(multiple testing correction: ', correction, ')', sep = '') );

  qq <- q[q[,1]<alpha,];
  if ( nrow ( qq ) > 0 ) {
    par(col='red',new=T);
    points(qq[nrow(qq):1,3],cutoff+1-nrow(qq):1,pch=16);
  }

  if ( src2 != "" ) {
    abline(v=mean(Mgood) - mean(M2good),col='blue');
  } else {
    if ( length ( classone ) == 0 && length ( classtwo ) == 0 ) {
      abline(v=controlmean,col='orange');
      abline(v=controlmean+controlsd,col='blue');
      abline(v=controlmean-controlsd,col='blue');
    } else {
      abline(v=mean(Mgood[,classone]) - mean(Mgood[,classtwo]),col='blue');
    }
  }

	dev.off();
  write.table ( signif(opvals[opvals[,1]<alpha,],3), sep="\t", quote=F, row.names = T, col.names = F, file = outputfile );
}


